Method for the determination of local similitude from seismic 3d measured data

ABSTRACT

The invention concerns a method for the determination of local similarity values for geological units in the subsurface from a seismic 3-D dataset, which consists of a multitude of traces, each of which being formed by a sequence of data points that carry amplitude values, or seismic attributes derived thereof, with the following steps: Assigning a weighted environment to a respective analysis point, whereby weight factors according to a weight function are assigned to the data points in the environment of the analysis point, and calculating similarity values in each weighted environment, whereby the weight function is included in the similarity measure; and especially a method for the determination of local dip-dependent similarity values with the steps: Assigning a weighted environment to a respective analysis point, whereby weight factors according to a weight function are assigned to the data points in the environment of the analysis point, calculating similarity values in each weighted environment, whereby the weight function is included in the similarity measure, and similarity values are determined for discrete spatial orientation, and the similarity values are represented at least as a function of a dip angle and a dip azimuth; and, in each weighted environment of an analysis point, determining the maximum similarity value, which is assigned, together with the corresponding dip angle and dip azimuth, to the respective analysis point.

BACKGROUND OF THE INVENTION

[0001] 1. Field of Invention

[0002] The invention concerns a method for the determination of localsimilarity of geological units in the subsurface from a seismic 3-Ddataset, which consists of a multitude of traces, each of which beingformed by a sequence of data points that carry amplitude values, orseismic attributes derived therefrom.

[0003] 2. Related Art of the Invention

[0004] Seismic reconnaissance methods are used worldwide in order toobtain additional knowledge on the distribution of geological structuresin the subsurface, in addition to the information from drilled wells.The information from seismic data often allows to renounce on furtherreconnaissance wells, or to reduce their number to a minimum.

[0005] For the seismic reconnaissance of the subsurface, sensors(geophones/hydrophones) are used. Being aligned in a line (2-Dseismics), they receive sound waves. These waves are emitted by aseismic source, e.g. by an explosive charge, vibrators, or airguns, andthey are partially reflected back to the surface by subsurface layers.There they are registered, and recorded in form of a time series. Thistime series represents the arriving seismic energy in form of amplitudefluctuations. It is stored digitally, and consists of regularly spaceddata samples, which are characterized by the respective time and by thecorresponding amplitude value. Such a time series is called a seismictrace. The recording line moves across the area to be examined, and thusthis a 2D seismic profile is recorded with this geometry.

[0006] The subsequent processing aims at a noise suppression, e.g. bystacking, or by dedicated application of filters. After stacking, whichsums up those reflection amplitudes that are assigned to the samesubsurface points, the term ‘poststack data’ is used. The results arevertical profiles, with a time or depth dependent representation ofamplitudes and traveltimes, as well as attributes derived from theamplitudes, which serve as a basis for further geologic evaluation. Thegeologic layers can be followed on a profile by the lateral alignment ofamplitudes.

[0007] If data are not recorded along a line but in an areal grid, athree-dimensional data volume results. In the case of a 3-D volume, anamplitude values is assigned to an arbitrary point in the subsurface,which may for example be described by Cartesian coordinates. Thevertical direction is measured in time (traveltime of sound).

[0008] In this case, a large amount of data is produced (severalGigabytes), which are stored and subjected to the processing, before theactual interpretation e.g. for the further reconnaissance of thesubsurface is possible. These processes require large computer resourcesand software, in order to process and correct the recorded signal. Theresult is given by poststack data, that represent a seismic volume inform of a 3-D dataset. The 3-D dataset represents physical andstructural characteristics of the subsurface in a seismic image.

[0009] From this dataset, arbitrary cuts may be extracted, e.g. verticalprofiles, and horizontal maps from different depths, which areinterpreted in the following by geophysicists and geologists. Since thisinterpretation of the achieved seismic images mainly comprises an opticcorrelation, attempts have been performed to automatize this subjectiveevaluation that depends on one, or more interpreters.

[0010] Accordingly, the U.S. Pat. Nos. 5,563,949 and 6,092,026 describea method, which highlights faults and zone of low coherency in athree-dimensional volume of stacked seismic data. This data volume issubdivided into a number of horizontal slices, and these slices areagain subdivided into a number of cells. The cells contain portions ofat least three seismic trace each, in a horizontal configuration thatallows a comparison in two pre-defined vertical planes, e.g. along theinline and crossline directions. As mathematical methods for determiningthe similarity, or coherency of the traces in these planes, thecross-correlation and the covariance are explicitly mentioned. Themaxima of the cross-correlations allow an estimation of the partial dipin the respective planes. These maxima can be determined for bothcorrelation planes, and be combined into a single coherency value by amathematical Operation. To each processed cell, the correspondingcoherency value is assigned, and a new seismic volume of the coherencyderived in this way, is thus created.

[0011] The calculation of trace similarity and dip by comparison of twotraces in inline direction, and in crossline direction, respectively, isvery sensitive to noise. Large time windows have to be chosen in orderto obtain stable results, thus reducing the resolution. The methodmeasures the local similarity/dissimilarity in cells, that possess afixed form and orientation with respect to the coordinates of theseismic data volume. A continuous, distance dependent correlation, andpossibly as well direction dependent correlation of the data can not beconsidered. The seismic data contained in the cell is entered into themeasurement of local similarity/dissimilarity with a constant weight.

[0012] A further method is known from the U.S. Pat. No. 6,141,622, whichmeasures the local similarity, or dissimilarity, respectively, of theseismic data in a three-dimensional seismic data volume. The measurementis carried out in a cell with vertical extension in time, and horizontalextension in the inline-crossline plane.

[0013] The traces used for the measurement are located either on oneline, i.e., on an inline or crossline, or they are located on two lines,i.e., on an inline and on a crossline, symmetrically to the intersectionof both lines. On these lines, 3, 5, or 7 neighboring traces areselected. The semblance, or the inverted semblance is measured in thecells. The measurement is performed exclusively along the inline, orcrossline included. For the two directions contained the x-shaped cell,two measurements are obtained, that are then summed together.

[0014] The method measures the local similarity/dissimilarity in cells,that possess a fixed form and orientation with respect to thecoordinates of the seismic data volume. A continuous, distance dependentcorrelation, and possibly as well direction dependent correlation of thedata can not be considered. The seismic data contained in the cell isentered into the measurement of local similarity/dissimilarity with aconstant weight. Dip is not considered.

[0015] The U.S. Pat. No. 6,138,075 describes a method that measures thelocal similarity, or dissimilarity of seismic traces in athree-dimensional seismic data volume. For the measurement, a cell isdefined in the vertical direction along the traces, and in thehorizontal plane. The cell contains a reference trace, and at least twoneighbor traces. For each neighbor trace, an individual value of thesimilarity to the reference trace is determined by applying mathematicalmeasures like semblance, or cross-correlation. The local dips of thedata volume are taken into account by vertically shifting the neighbortraces in a pre-defined range, and selecting the maximum similarityvalue. The neighbor trace with the largest similarity value is declaredthe target trace. Then all similarity values that have been derived upto this stage, are regarded as preliminary, and are deleted. The finalsimilarity measurement takes place exclusively between

[0016] The reference trace and the target trace, and it may well use adifferent similarity measure, e.g. the ‘Manhattan distance’.

[0017] The determination of similarity values by comparison of twotraces, respectively, is very sensitive to noise. The limitation of thecomparison to neighbor traces is not suited to detect gradual changes ofthe seismic behavior, which extend over distances of the order of manytrace intervals, and where the change in the range of one trace intervalis below the noise level of the data.

[0018] Another method is known from the U.S. Pat. No. 6,151,555, or WO00/54207, respectively, which measures the local dissimilarity ofseismic traces in a three-dimensional data volume. Tow alternativemeasures of the dissimilarity are given, following the concept of thestatistical variance. For the measurement, a cell is defined in thevertical direction along the traces, and in the horizontal plane. Forthe horizontal plane. quadratic cells with 3×3 or 5×5 traces, as well asan x-shaped cell with 5+5 traces in both orthogonal directions, areoffered for selection. The vertical extension of the cell is describedby the number of contained samples along each single trace; thissimultaneously represents the number of horizontal data slices. Atriangular weight function assigns a weight to each horizontal dataslice. Within the horizontal plane, however, no weighting takes place.Hence, each horizontal position contributes in the same way to themeasure of dissimilarity. In each horizontal data slice, the mean, thesum of the squared deviations from the mean, and the sum of the squaresis derived from data values contained in the slice. Two measures ofdissimilarity can be calculated from these quantities as follows:

[0019] 1. The sums of the squared deviations are multiplied with thevertical weight function, and summed for all horizontal data slices; thesums of the squares are treated equivalently, and finally the former aredivided by the latter.

[0020] 2. Alternatively, an individual quotient is calculated for eachhorizontal data slice, by dividing the sum of the squared deviations bythe sum of the squares; these quotients are multiplied with the verticalweight function, and finally summed for all horizontal data slices.

[0021] The method measures the local dissimilarity in cells, whichpossess a fixed form and orientation with respect to the coordinates ofthe seismic data volume. A continuous, distance dependent correlation,and possibly as well direction dependent correlation of the data can notbe considered. The seismic data contained in the cell is entered intothe measurement of local similarity/dissimilarity with a constantweight. A consideration of dip does not take place.

[0022] Furthermore, a method is known from the DE 199 33 717 C1 whichrelates to the similarity analysis of data points belonging to apre-defined three-dimensional environment of the analysis point justconsidered in a three-dimensional seismic dataset. The three-dimensionalseismic dataset consists of a multitude of traces, each of which iscomposed as a sequence of data points to which amplitude values areassigned, or seismic attributes derived from amplitude values. Thismethod, however, calculates the similarity of local sections of seismicdata from the measurement dataset to a reference section correspondingto a pre-defined location and depth, and assigns this similarity valueas an attribute to the central value of the corresponding local section.

[0023] A method is known from the WO 97/13166, which measure the localdissimilarity of seismic traces in a three-dimensional seismic datavolume with a semblance method. For the measurement, a cell is definedin the vertical direction along the traces, and in the horizontal plane.In the horizontal plane, the cell has an elliptical, or rectangularshape, and it is centered at the analysis point. The main axis of thecell may exhibit an arbitrary angle to the horizontal co-ordinate axisof the seismic data. However, its orientation with respect to the datagrid is then rigidly fixed by this angle. The measurement of thesimilarity may thus include one direction of preference at most, whichis defined by the direction of the main axis. The vertical extension ofthe cell is described by the number of contained samples along eachsingle trace; this simultaneously represents the number of horizontaldata slices.

[0024] The semblance is calculated with in the cell for complex traces,i.e., for each real trace, a Hilbert transformed imaginary trace has tobe derived. An estimation of the coherency follows from the averaging ofthe semblance over the vertical extension of the cell. In eachhorizontal data slice, the sum of the real values, and the sum of theHilbert transformed values are calculated. Each of these sums issubsequently squared. These squared amplitude sums of all horizontaldata slices are summed in order to form the numerator of the semblance.The denominator contains the sum of the squares of all real, and Hilberttransformed single values that are contained in the cell. The semblancein calculated as a dip dependent quantity, with dip being described bytwo parameters. These parameters are the apparent dips in the directionsof the orthogonal horizontal data axes. They can be transformed into theparameter pair of dip and azimuth angles. For the dip dependentsemblance calculation, the vertical position of the cell in thethree-dimensional data volume is changed, whereas its horizontal shapeand centering at the analysis point is preserved. The cell is shearedaccording to the apparent dips, i.e., the partial sums for the semblancecalculation are derived in correspondingly dipping data slices. The dipsto be considered are limited by a maximum dip angle. Three alternativeschemes for the discretization of the solid angle are given in a polarrepresentation of the dip and azimuth angle. The discretization may beperformed with a quadratic, a triangular, or a radial grid. The radialgrid is hereof regarded as disadvantageous.

[0025] Contrasting to the other schemes, it involves a highlynon-uniform discretization of the solid angle.

[0026] The presented method does not use weighting in the horizontalplane, which could take into account a decrease of the lithologic andstructural correlation of geologic bodies with distance. A temporalweighting is not applied either. The horizontal cell may possess adirection of preference, the orientation of which remains unchangedduring the whole process of the similarity or coherency determination inthe three-dimensional data volume. Hence, during the dip measurement,this rigidly oriented cell is not oriented according to the dip azimuth,i.e. to the direction of structural strike that varies more or less inthe considered volume. If the directions of geologic strike, which areassumed during the dip dependent similarity measurements, vary withrespect to the direction of preference of the cell, this may bedisadvantageous, since the direction of largest continuity in geologicalbodies often has a relation to the direction of geological strike.Consequently, a possibly existing direction of preference of the cellmay be oriented partly in parallel, partly perpendicular, and frequentlywith other angles to the directions of geologic strike of the considereddips. However, here it is possible that the measured similarity valuesrender a distorted view of the exising dip dependence in the data, sincethe measurement renders systematically increased, or decreasedsimilarity for some azimuth ranges due to the orientation of theanalysis cell.

[0027] Moreover it is disadvantageous, that the Hilbert transformedimaginary trace must be calculated for each real trace within thecalculation of similarity with consideration of dip and dip azimuth.This additionally requires significant calculation capacities. Thismethod known from the WO 97/13166 is also revealed with supplementaryexamples in the publication by Marfurt, K. J., et al.: 3-D seismicattributes using a semblance-based coherency algorithm.—In: Geophysics,1998, volume 63, page 1150-1165.

[0028] All previously mentioned, known methods perform the measurementof a similarity value in cells within a three-dimensional volume ofseismic poststack data. The cells possess a fixed form and size in thespatial horizontal plane. The partial sequences of seismic traces thatare contained in a cell, contribute with a horizontally constant weightto the derivation of individual similarity values for single pairs oftraces, or of a general similarity value for the whole cell. Hence, thetrace weighting jumps from a constant positive level inside the cell tothe value zero outside the cell. This means, that all methods cannot thetake into account the degree of spatial continuity of geological bodies,as it is derived e.g. by variographic analysis, since neighborhoodrelations are not weighted. Variography calculation and modeling is aprerequisite of geostatistics. It is used to describe quantitatively thedistance dependent, and, if appropriate as well direction dependentspatial relations of neighboring points.

SUMMARY OF THE INVENTION

[0029] The purpose of the invention is to provide a processing methodfor seismic 3-D data for the determination of the local similarity ofgeological units in the subsurface, where in this method the spatialcontinuity is considered with a distance dependence, and, ifappropriate, as well with a direction dependence.

[0030] This task is solved by a method according to claim 1.

[0031] Hence, the invention comprises a method for the automaticdetection and evaluation of structural, facial and lithologic units withdiscontinuities and transitions of different spatial magnitudes inseismic data. It is designed for the processing of a multitude ofseismic traces from a three-dimensional data volume. Thethree-dimensional data volume can be arranged in a way, that thelocation of a single trace is characterized by two coordinates (x,y) ina horizontal plane, whereas the third, vertical coordinate axis isdirected along the vertically arranged seismic traces. The thirdcoordinate (z) can describe a time or depth on a seismic trace. Fordiscrete coordinate values (x_(i), y_(j),z_(k)), the three-dimensionalvolume contains single data values, i.e., samples s(x_(i),y_(j),z_(k))of the seismic amplitude, or of a seismic attribute derived thereof.

[0032] By the weighted environment of each analysis point, a spatiallyweighted determination of the neighborhood relation is provided for thedetection of facial, lithologic and structural units, and theirtransitions, respectively. The dimension and weighting of theenvironment that is used for the determination of the neighborhoodrelation can be selected in analogy to the spatial relations determinedby variographic analysis. Alternatively, a weight function can beextracted from the scientific literature with dedication to certaingeologic problems. Weighted environments of any form and size can beused for similarity measurements, since the weight function is usuallydefined on the whole data volume. For the similarity evaluation, theweighted environment can be combined with an arbitrary similaritymeasure, by incorporating the weight function into the selectedsimilarity measure. Common similarity measures for seismic data are e.g.the semblance, the correlation, or the variance.

[0033] This similarity analysis is preferably carried out withconsideration of the local dips, by calculating dip dependent similaritymeasures according to claim 2.

[0034] The reflections of a three-dimensional seismic data volumerepresent the layer planes, and tectonic boundary planes of thesubsurface, which often do not extend horizontally, but with dips ofvarying strengths and directions.

[0035] Correspondingly, dipping reflections occur in the seismic datavolume as well.

[0036] For such reflections with varying dips, comparable similarityvalues must be determined. For this purpose, the similarity measurementis carried out along the local dip of the seismic data. For this dipdependent similarity measurement, the weighted environment is shearedvertically, such that the previously horizontal planes of thisenvironment are inclined according to the local dip.

[0037] However, the local dip of the seismic data at an analysis pointis unknown at the beginning of the similarity measurements, and must beestimated. Here it is assumed that the similarity measurements atdifferent dips render a maximum similarity value in that case, where thedip of the measurement, and the local dip in the data volume coincide.

[0038] In order to find this dip-dependent maximum with sufficientaccuracy, the range of possible dips is subdivided according to apreferably uniform discretization of the solid angle.

[0039] After the detection of the dip dependent coherency maximum, thiscoherency maximum and the corresponding dip and azimuth angles arestored as the results of the similarity measurements, and assigned tothe respective analysis point. They are then available for furtherdigital evaluation, e.g. by pattern recognition. Furthermore, they aregraphically displayed to the interpreter on the monitor screen, alongcuts through the three-dimensional volume, or along the evaluated timeslices, or horizons, respectively. Clearly, these displays can beprinted as well.

[0040] If the weight function is centered at the analysis point, andexhibits decreasing values away from the analysis point, then theinfluence of a data value on the similarity measure will decrease withincreasing distance from the analysis point, as has to be commonlyassumed.

[0041] If the weighting is embodied by window functions, especially bytriangular, Hamming, Hanning or Daniell functions, then an environmentthat is limited by the window function is chosen for the similarityevaluation.

[0042] By using a weight function that decays with distance (d) from theanalysis point up to infinity, especially by using d^(−c), ^(−cd),e^(−cdd) with c=arbitrary constant, and by setting the weight to zerobeyond a maximum distance (d_(max)) from the analysis point, thecalculation effort for functions that decay for an infinite distancerange, is limited as well. At the maximum distance from the analysispoint, the transition to the weight factor zero can be smoothed, e.g. byan additional linear transition function, in order to reduce edgeeffects in the calculation of similarity values.

[0043] A finite environment for the similarity evaluation is alsoprovided by an arbitrary positive weighting, where a window functioncauses a tapering to the value zero in an edge area. Here the weightingis formed as the product of two functions, where one of the functionsdoes not decay, or does not fully decay with distance from the analysispoint, and the other function represents a window function, whichpossesses the value 1 in the central area of the environment, decreasesfrom 1 to 0 in an edge area of the environment, and exhibits the value 0outside the environment.

[0044] By forming a weighted environment with a rotational symmetryaround the analysis point at least in the horizontal plane, adirectionally neutral environment is defined, which does not exhibit adirection of preference in the similarity calculations.

[0045] A distance dependent weighting can be equipped with a directionof preference by a direction dependent scaling of the distance measure,preferably by an elliptical scaling. By rotating a weighted environmentthat does not exhibit a rotational symmetry around the analysis point,according to the respective dip azimuth for the determination ofsimilarity values, and especially by orienting an existing direction ofpreference that is due to the non-symmetry, in a fixed geometricalrelation with respect to the dip azimuth, it can be taken into accountthat in geological bodies the direction of largest continuity is oftenrelated to the direction of geological strike, or to the dip azimuth,respectively. By rotating the environment that does not exhibit arotational symmetry, together with the respective discrete dip azimuthsthat are considered in the measurement of similarity, a distortion ofthe similarity values can just be avoided for typical geologicalstructures.

[0046] If the analysis interval between neighboring analysis points inthe 3-D dataset is selected as a multiple, and especially as an integermultiple of the data point interval, preferably between one and severaltens, especially preferred between three and ten such data pointintervals, the computation time can be significantly reduced, since asimilarity measurement must not be performed for every sample position.The selection of the analysis points is defined by lateral and verticalanalysis intervals. These analysis intervals are preferably integermultiples of the data point intervals, i.e., the sample intervals.However, rational multiples are as well possible, requiring anassignment of interpolated data values to positions that are definedbetween the measured data points. After calculating a similarity valueand a dip with corresponding azimuth at an analysis point, the usedweighted environment is centered at the analysis point that is nextaccording to the pre-defined analysis interval. For environments with alarge spatial extension, large analysis intervals between neighboringanalysis points can generally be used in order to reduce the computationeffort, without deteriorating the resolution.

[0047] If the dip and azimuth angles are discretized at constant radialand angle intervals in a radial grid, the solid angle can be classifiedwith advantageous simplicity in the statistical evaluation, and in colordisplays of strike directions. However, it is disadvantageous that theconstant angle intervals imply a very non-uniform discretization of thesolid angle. Near the grid origin (horizontal layering), the analysispoints are narrowly spaced, which means an unnecessary load on thecomputing capacity.

[0048] Hence, it is proposed as a further development, that, for thecalculation of similarity values at discrete spatial orientations, thediscrete dip angle (θ_(i)) is calculated as follows, starting from apre-defined maximum dip angle (θ_(max)): at dip angles θ_(i) withθ_(max)/2<θ_(i)≦θ_(max), a similarity value is calculated for everydiscretized azimuth, at dip angles θ_(i) with θ_(max)/4<θ_(i)≦θ_(max)/2,the calculation is performed for every second azimuth, at dip anglesθ_(i) with θ_(max)/8<θ_(i)≦θ_(max)/4 for every forth azimuth, at dipangles θ_(i) with θ_(max)/16<θ_(i)≦θmax/8 for every eighth azimuth,etc., with a preferred termination at the fourth subdivision. Thus, thedisadvantage of the non-uniform discretization of the solid angle in theradial grid, and of the corresponding computation effort for manysimilarity determinations near the grid origin, is compensated for.However, the discretization of the dip and azimuth angles in constantsteps, which is advantageous for further evaluations, is retained. Inthe dip estimation as described before, the accuracy of thedetermination of the azimuth angle can be refined during an optionalfurther stage, if for a dip estimate in a coarsened azimuth interval,similarity values are calculated at the original angle steps on bothsides of the first estimate, as long as an increase of the similarityvalues is observed. The first dip estimate is replaced by the dip andazimuth angles of the general maximum, if the maximum of the similarityvalues from this second stage exceeds the maximum from the first stage.

[0049] When computing the similarity along arbitrarily oriented planesas the square of the representative attribute (amplitude), divided bythe average of the squared attribute (amplitude), the representativeamplitude replaces the mean amplitude of the semblance measure which iscommon in seismics. The representative amplitude is computed, e.g. by aneuronal network.

[0050] Preferably, the representative attribute is the weighted medianof the attributes to be analyzed. The application of the weighted medianof the amplitude (the attribute) implies a sorting out of single strongspikes and noise amplitude, which in case of the semblance can have astrong influence on the result. This sorting-out improves the analysisresult.

[0051] According to claim 14, the weighted median coherency is modifiedfor the application with weighted environments.

[0052] The weighted coherency on the basis of the weighted semblance, asgiven in claim 15, does as well take into account the weight factors ofthe similarity determination. In the following, an example of anembodiment is described in details, with reference to the encloseddrawings.

[0053] It is illustrated:

[0054]FIG. 1 a discretization of the solid angle by a modified radialgrid,

[0055]FIG. 2 for a time slice of a seismic in a) the obtained maximumlocal similarity value, in b) the corresponding dip angle and in c) thecorresponding dip azimuth.

[0056] First, the basis for calculation and data evaluation according tothe invention is described in the following.

[0057] Weight Function

[0058] Around each analysis position (x_(I),y_(J),y_(K)) that isselected for the calculation of a similarity value, a three-dimensionalenvironment is defined by a weight function

g _(IJK)(x _(i) ,y _(j) ,z _(k))=g(x _(i) −x _(I) ,y _(j) −y _(J) ,z_(k) −z _(K)).   (1)

[0059] This weight function contains arbitrary values g(x,y,z)>0 in athree-dimensional environment of arbitrary shape and arbitrary extent,around the co-ordinate origin (x,y,z)=(0,0,0), and it contains valuesg(x,y,z)=0 outside this environment. It contains the weights, which areimposed on the seismic data values located in the environment, when theyare entered into the calculation of a similarity value at the selectedanalysis location (x_(I),y_(J),z_(K)).

[0060] This function can e.g. exhibit decaying weights at the edge ofthe environment, and thus reduce edge effects in the calculation ofsimilarity values. Just as well, a stronger weighting of the valueslocated closer to the origin, can improve the result, since this favorsthe direct neighbor values in the calculation of the similarity value.Generally, this weighting can also be derived from variograms, whichdescribe the distance dependent correlation of parameters in theinvestigation area.

[0061] Special weight functions can be composed from independent partialfunction for single co-ordinates, e.g. in the form

g(x,y,z)=p(x,y)ƒ(z),   (2a)

or

g(x,y,z)=l(x)h(y)ƒ(z),   (2b)

[0062] where ƒ(z) describes an arbitrary function for the extraction oftime or depth windows from seismic traces. Here, for example, atriangular temporal weighting can be used in combination with ahorizontal weighting p(x,y), or l(x)h(y), according to the equations(2a) and (2b).

[0063] The horizontal weight function p(x,y) in equation (2a) can beselected e.g. in analogy to the spatial relation derived by variographicanalysis. From distance dependent, and, if appropriate, directiondependent values of the spatial correlation, value tables can bederived, which are subsequently interpolated, or approximated byanalytical functions. Spatial directions of preference can beemphasized, or reduced along the horizontal coordinate axes by scalingfactors a,b>0 according to

g(x,y,z)=p _(a,b)(x,y)ƒ(z)=p(ax,by)ƒ(z),   (3a)

[0064] and be rotated in the horizontal plane by an angle α according to

g(x,y,z)=p _(a,b,α)(x,y)ƒ(z)=p(a[x cos α−y sin α], b[x sin α+y cosα])ƒ(z)   (3b)

[0065] From this, special weightings with elliptical and circulargeometry follow with the elliptic equation

r _(a,b,α)(x,y)=(a ² [x cos α−y sin α]² ,b ² [x sin α+y cosα]²)^(0.5)=const   (4a)

as

g(x,y,z)=q(r _(a,b,α)(x,y))ƒ(z).   (4b)

[0066] The circular geometry follows with a=b.

[0067] The function q(d) in equation (4b) describes a distance dependentweighting. Here, weightings are preferred which decay with distance.Every window function that is common in seismic data processing may beused, e.g. the triangular, Hamming-, Hanning-, Daniell functions, orother weight functions. Beyond a maximum distance (d>d_(max)), thewindow functions usually adopt the value zero.

[0068] In case functions that decay for values up to infinity, thecomputation effort can be limited as well by using these functions onlyup to a maximum distance d_(max), and setting the function value to zeroat larger distances (d>d_(max)). At d_(max), a jump of the weightfunction q(d) is thus created, which can be compensated for by localsmoothing. Suitable functions that decay with distance d up to infinityare d^(−c), e^(−cd), e^(−cdd), and other functions.

[0069] The functions q(d) listed here, can also be used for a weightingaccording to equation (2b) with l(x)=q(|x|) or h(y)=q(|y|).

[0070] As long as no other information favors the selection of a specialweight function, the common practice has shown that the similaritycalculation with the Gaussian bell shaped exponential curve as ahorizontal weight function renders good results. The environment of ananalysis location in equation (1) is defined by a special weightfunction according to equation (4a,b). In this way, the distancedependent weighting has the form $\begin{matrix}{{q(d)} = \begin{Bmatrix}^{- d^{2}} & {for} & {0 \leq d \leq 1.5} \\{0.5797 - {0.3162\quad d}} & {for} & {1.5 < d \leq 1.8333} \\0 & {for} & {1.8333 < d}\end{Bmatrix}} & ( {5a} )\end{matrix}$

[0071] with the scaled distance

d=r _(a,a,0)(x,y)=a(x ² +y ²)^(0.5).   (5b)

[0072] For a similarity determination with high resolution, a distancescaling according to a⁻¹=2Δu has proved to be efficient, with Δu beingthe average interval of the horizontal discretization in the seismicvolume. For the resolution of facies changes, the reciprocal scalingfactor a⁻¹ should be chosen in the order of characteristic lengths ofthe geologic bodies.

[0073] Similarity Measure

[0074] In connection with a weighted environment, most of the commonsimilarity measures have to be modified. With help of suchmodifications, however, all known similarity measures can be applied aswell in horizontally and vertically weighted environments. A significantmodification is the following, that each summed value is weighted in thesums over amplitudes, squared amplitudes, or other functions of theamplitude, and that the normalization of such sums does not incorporatethe number of summed values, but the sum of the corresponding weightfactors.

[0075] For two special similarity measures, this modification isdescribed in the following:

[0076] The weighted environment g_(IJK)(x_(i),y_(j),z_(k)) of ananalysis location (x_(I),y_(J),z_(K)) is defined in equation (1). Asimilarity measurement is performed within this environment along ahorizontal data slice at an arbitrary time or depth z_(k). In analogy tothe semblance that is common in similarity measurements, the followingweighted semblance can be used: $\begin{matrix}\begin{matrix}{{S_{\overset{\_}{k}}( {x_{I},y_{J},z_{K}} )} =} \\\frac{\{ {\sum\limits_{i,j}{g_{IJK}( {x_{i},y_{j},z_{\overset{\_}{k}}} )}} \}^{- 1}\{ {\sum\limits_{i,j}{{g_{IJK}( {x_{i},y_{j},z_{\overset{\_}{k}}} )}{s( {x_{i},y_{j},z_{\overset{\_}{k}}} )}}} \}^{2}}{\sum\limits_{i,j}{{g_{IJK}( {x_{i},y_{j},z_{\overset{\_}{k}}} )}\{ {s( {x_{i},{y_{j}z_{\overset{\_}{k}}}} )} \}^{2}}}\end{matrix} & ( {6a} )\end{matrix}$

[0077] By summation over the entire time or depth range of the weightedenvironment, the similarity measure of the weighted coherency results:$\begin{matrix}\begin{matrix}{{S( {x_{I},y_{J},z_{K}} )} =} \\\frac{\sum\limits_{k}\lbrack {\{ {\sum\limits_{i,j}{g_{IJK}( {x_{i},y_{j},z_{k}} )}} \}^{- 1}\{ {\sum\limits_{i,j}{{g_{IJK}( {x_{i},y_{j},z_{k}} )}{s( {x_{i},y_{j},z_{k}} )}}} \}^{2}} \rbrack}{\sum\limits_{k}\lbrack {\sum\limits_{i,j}{{g_{IJK}( {x_{i},y_{j},z_{k}} )}\{ {s( {x_{i},{y_{j}z_{k}}} )} \}^{2}}} \rbrack}\end{matrix} & ( {6b} )\end{matrix}$

[0078] With the normalization in the nominator as in equation (6b),using the sum of the weights, the similarity measure is in many casesmore uniformly influenced by the amplitude distribution of the wholeenvironment, than with the normalization in the denominator as describedin the following: $\begin{matrix}\begin{matrix}{{\overset{\_}{S}( {x_{I},y_{J},z_{K}} )} =} \\\frac{\sum\limits_{k}\{ {\sum\limits_{i,j}{{g_{IJK}( {x_{i},y_{j},z_{k}} )}{s( {x_{i},y_{j},z_{k}} )}}} \}^{2}}{\sum\limits_{k}\lbrack {\{ {\sum\limits_{i,j}{g_{IJK}( {x_{i},y_{j},z_{k}} )}} \} {\sum\limits_{i,j}{{g_{IJK}( {x_{i},y_{j},z_{k}} )}\{ {s( {x_{i},y_{j},z_{k}} )} \}^{2}}}} \rbrack}\end{matrix} & ( {6c} )\end{matrix}$

[0079] The weighted coherency in the formulation (6b), or in themodification (6c), contains two-step summations: The inner summationsconcern the positions in a horizontal plane, the outer summations, onthe contrary, concern the horizontal planes that are present in theweighted environment. This is an orientation sensitive measurement ofsimilarity, which measures similarity in the horizontal direction. Anorientation insensitive similarity measurement generally comprises onlyone-step summations over all data points that are contained in theweighted environment. In this way, the weighted coherency is obtained asan orientation insensitive similarity if in equation (6a), each of thesummations over the horizontal indices i,j is extended by the summationover the vertical index {overscore (k)}.

[0080] The seismic amplitude values in an environment do on one handreflect the geology of the subsurface, but on the other hand they arepartly as well the result of a disadvantageous signal-to-noise ratio.Hence, this invention prefers a similarity measurement, which improvesthe signal-to-noise ratio in the similarity determination.

[0081] This concept is based on a generalization of the semblance, andcomprises the weighted median similarity, and the median coherency, asdescribed in the following. The application of these similarity measuresis performed in analogy to the similarity measures previously described.

[0082] Again, a horizontal data slice at an arbitrary time or depthz_({overscore (k)}) with a total of N amplitude valuess(x_(i),y_(j),z_({overscore (k)})), and with N corresponding weightfactors g_(IJK)(x_(i),y_(j),z_({overscore (k)})) is considered. Theamplitude values are sorted in increasing order, and newly indexed:

s₁≦s₂≦s₃≦ . . . ≦s_(N−1)≦s_(N).

[0083] This sorting and indexing are as well transferred to the valuesof the weight function g_(IJK), that belong to the seismic amplitudes:

g, g₂, g₃, . . . g_(N−1), g_(N)

[0084] From this, partial sums of the weight factors are derived:

G ₁ ≦G ₂ ≦G ₃ ≦ . . . ≦G _(N−1) ≦G _(N), where G _(n) =Σg _(i)(i=1 ton).

[0085] Now the weighted median can be derived as $\begin{matrix}\begin{matrix}{{m_{IJK}( z_{k} )} =} \\{\begin{Bmatrix}{s_{n},} & {{{{if}\quad G_{n - 1}} < \frac{G_{N}}{2}},{G_{n} > \frac{G_{N}}{2}}} \\{\frac{s_{n} + s_{m}}{2},} & {{{{if}\quad G_{n - 1}} < \frac{G_{N}}{2}},{G_{n} = {\frac{G_{N}}{2} = G_{m - 1}}},{G_{m} > \frac{G_{N}}{2}}} \\\quad & {{{with}\quad n} < m}\end{Bmatrix}.}\end{matrix} & ( {7a} )\end{matrix}$

[0086] In order to perform a similarity measurement within the weightedenvironment along a horizontal data slice at an arbitrary time or depthz_(k), the weighted median similarity is formulated as follows with thepreviously defined weighted median: $\begin{matrix}{{M_{\overset{\_}{k}}( {x_{I},y_{J},z_{K}} )} = \frac{\{ {\sum\limits_{i,j}{g_{IJK}( {x_{i},y_{j},z_{\overset{\_}{k}}} )}} \} \{ {m_{IJK}( z_{\overset{\_}{k}} )} \}^{2}}{\sum\limits_{i,j}{{g_{IJK}( {x_{i},y_{j},z_{\overset{\_}{k}}} )}\{ {s( {x_{i},y_{j},z_{\overset{\_}{k}}} )} \}^{2}}}} & ( {7b} )\end{matrix}$

[0087] From this, the weighted median coherency follows as a similaritymeasure by summing over the whole time or depth range of the weightedenvironment: $\begin{matrix}{{M( {x_{I},y_{J},z_{K}} )} = \frac{\sum\limits_{k}\lbrack {\{ {\sum\limits_{i,j}{g_{IJK}( {x_{i},y_{j},z_{k}} )}} \} \{ {m_{IJK}( z_{k} )} \}^{2}} \rbrack}{\sum\limits_{k}\lbrack {\sum\limits_{i,j}{{g_{IJK}( {x_{i},y_{j},z_{k}} )}\{ {s( {x_{i},{y_{j}z_{k}}} )} \}^{2}}} \rbrack}} & ( {7c} )\end{matrix}$

[0088] Like the weighted coherency in the formulations (6b) and (6c),the weighted median coherency in equation (7c) represents an orientationsensitive similarity measurement, which measures the similarity withinthe weighted environment in horizontal direction. In analogy to theremarks concerning the weighted coherency, the weighted median coherencymay as well be formulated as an orientation insensitive similaritymeasure if in equation (7b), each of the summations over the horizontalindices i,j is extended by the summation over the vertical index{overscore (k)}, and thus covers all positions within the weightedenvironment. Moreover, the weighted median m_(IJK) in equation (7b) mustnot be calculated for a horizontal data slice at an arbitrary time ordepth z_({overscore (k)}), but for the entire weighted environment. Thisis accomplished in analogy to equation (7a) by sorting all amplitudevalues and corresponding weight factors that are contained in theweighted environment.

[0089] Horizontal Orientation, and Dip of the Weighted Environment

[0090] The layer planes and tectonic boundary planes of the subsurfaceare often not oriented horizontally, but possess dips of varyingstrengths and directions. Accordingly, the reflections in athree-dimensional seismic data volume possess different dips as well.Separate similarity values, being calculated in the seismic data volumefor local dips, are comparable, if the similarity measurement wascarried out with constant orientation with respect to the dip in allcases. Hence, for dip dependent similarity measurements in a weightedenvironment, the weighted environment is

[0091] 1. rotated according to the dip direction in the horizontalplane, i.e., is rotated for the respective similarity measurement, and

[0092] 2. inclined in space according to the dip.

[0093] For the similarity measurement in case of horizontal layering, orhorizontal orientation of the reflections, respectively, the weightedenvironment is described in general by equation (1), and in specialembodiments by the equations (2)-(5). This weighted environment isreferred to in the following as the basic environment, in which nosecondary rotations, or dips have been introduced. The basic environmentmay generally exhibit a horizontal direction of preference. Thisdirection of preference must retain a constant angle to the respectivedip direction in similarity measurements with different dip directions.Such an orientation of the basic environment in the horizontal plane isnot required only in cases, where no unique direction of preference ispresent, i.e., in cases of predominant of full horizontal symmetry.

[0094] The equations (1)-(5) define the basic environment with aninitial horizontal orientation. In case of dip dependent similaritymeasurements, the orientation of a weighted environment is tied to thedirection of dip. This is fixed by assigning a certain horizontal dipdirection, or dip azimuth φ_(o), respectively, to the initial horizontalorientation of the basic environment.

[0095] This dip direction is commonly selected in parallel, orperpendicular to the horizontal direction of preference; however, it canas well have an arbitrary angle to it. For any other dip directionΦ=φ_(o)+φ, the weighted environment is rotated by the angle φ, andsheared from the horizontal orientation according to the dip angle θ.Thus, a modified form of the weighted environment from equation (1) isobtained as

g _(IJK,φ,θ)(x _(i) ,y _(j) ,z _(k))=g _(φ,θ)(x _(i) −x _(I) ,y _(j) −y_(J) ,z _(k) −z _(K)).   (8a)

[0096] In absence of a direction of preference, a shearing according tothe dip angle θ in the dip direction Φ_(g)=φ_(o)+φ is sufficient:

g″ _(IJK,φ,θ)(x _(i) ,y _(j) ,z _(k))=g″ _(φ,θ)(x _(i) −x _(I) ,y _(j)−y _(J) ,z _(k) −z _(K)).   (8b)

[0097] In case of a dip in the direction of the basic dip azimuth φ_(o),only a shearing of the weight function according to the dip angle θresults. The weight function g_(φ,θ) of equation (8a) thus receives theform

g _(0,θ)(x,y,z)=g _(θ)(x,y,z)=g(x,y,z+[x sin φ_(o) +y cos φ_(o)] tan θ).  (9a)

[0098] In case of an azimuth Φ_(g)=φ_(o)+φ with φ>0°, the weightfunction with a direction of preference is additionally rotated in thehorizontal plane:

g _(φ,θ)(x,y,z)=g _(θ)(x cos φ−y sin φ, x sin φ+y cos φ, z).   (9b)

[0099] In absence of a direction of preference, a rotation of the weightfunction as in equation (9b) is not required. A shearing in the dipdirection Φ_(g) is sufficient:

g″ _(φ,θ)(x _(i) ,y _(j) ,z _(k))=g(x,y,z+[x sin Φ_(g) +y cos Φ_(g)] tanθ).   (9c)

[0100] The weight function g in the equations (9a,c) is definedaccording to the equations (1)-(5).

[0101] Similarity Measure in Case of Horizontal Orientation and Dip ofthe Weighted Environment

[0102] In an environment with directional orientation in the horizontalplane, and with dip, the similarity measurements as well take place indipping planes, that possess a dip angle θ and a dip azimuth φ. Thesemeasurements are performed in analogy to the basic environment withoutdip:

[0103] In the basic environment without dip g_(IJK)(x_(i),y_(j),z_(k))for an analysis position (x_(I),y_(J),z_(K)), individual parameters ofthe similarity measurement are determined in single horizontal planes atthe vertical co-ordinate z_({overscore (k)}). For special similaritymeasurements, this is demonstrated in the equations (6a), (7a), (7b) atthe parameters S_({overscore (k)}),m_(IJK)(z_({overscore (k)}))M_({overscore (k)}), where exclusively theseismic amplitudes s(x_(i),y_(j),z_({overscore (k)})) and correspondingweight factors g_(IJK)(x_(i),y_(j),z_({overscore (k)})) are used.

[0104] In a weighted environment with dip θ and dip azimuth Φ_(g)=0°,the individual dipping planes are characterized by their location indepth, or time, respectively:

z _({overscore (k)},φ,θ)(x _(i) ,y _(j))=z _({overscore (k)})−tan θ(y_(j) −y _(J))

[0105] For a dip azimuth Φ_(g)≠0°, the planes possess the depth location

z _(k,φ,θ)(x _(i) ,y _(j))=z _({overscore (k)})−tan θ└(x _(i) −x_(I))sin φ+(y _(j) −y _(J))cos φ┘

[0106] Here, z_(k) characterizes the intersection (x_(I),y_(J),z_(k)) ofthe individual plane with a vertical axis through the analysis position(x_(I),y_(J),z_(K)). The depth location z_(k,φ,θ) of the dipping planeusually does not coincide with the points of the vertical discretizationof the seismic volume. Hence, from the data samples s(x_(i),y_(j),z_(k))at the horizontal discretization point (x_(i),y_(j)), a seismicamplitude {tilde over (s)}(x_(i),y_(j),z_(k,φ,θ)) has to be interpolatedat the depth location z_(k,φ,θ)(x_(i),y_(j)).

[0107] The similarity measurements in horizontal planes can betransferred to dipping planes by two changes:

[0108] The discretized values s(x_(i),y_(j),z_(k)) of the seismicamplitude are replaced by values {tilde over(s)}(x_(i),y_(j),z_(k,φ,θ)), that are interpolated, if required.

[0109] A weight function g_(IJK)(x_(i),y_(j),z_(k)) with a direction ofpreference is replaced by g_(IJK,φ,θ)(x_(i),y_(j),z_(k)) according tothe equations (8a), (9b). A weight function g_(IJK)(x_(i),y_(j),z_(k))without a direction of preference, however, is replaced byg″_(IJK,φ,θ)(x_(i),y_(j),z_(k)) according to the equations (8b), (9c).

[0110] These changes are illustrated for the example of a specialsimilarity measure. In case of dip and rotation, the similarity measurethat is formulated in equation (6b) for the basic weighted environment,adopts the form $\begin{matrix}\begin{matrix}{{S_{\varphi,\vartheta}( {x_{I},y_{J},z_{K}} )} =} \\\frac{\sum\limits_{k}\lbrack {\{ {\sum\limits_{i,j}{g_{{IJK},\varphi,\vartheta}( {x_{i},y_{j},z_{k}} )}} \}^{- 1}\{ {\sum\limits_{i,j}{{g_{{IJK},\varphi,\vartheta}( {x_{i},y_{j},z_{k}} )}{\overset{\sim}{S}( {x_{i},y_{j},z_{k,\varphi,\vartheta}} )}}} \}^{2}} \rbrack}{\sum\limits_{k}\lbrack {\sum\limits_{i,j}{{g_{{IJK},\varphi,\vartheta}( {x_{i},y_{j},z_{k}} )}\{ {\overset{\sim}{S}( {x_{i},y_{j},z_{k,\varphi,\vartheta}} )} \}^{2}}} \rbrack}\end{matrix} & (10)\end{matrix}$

[0111] Determination of Local Dip

[0112] In the above, methods for the similarity measurements have beenpresented, which allow to consider local dip in the seismic data volume.Yet, the local dip is unknown. It is assumed, however, that at ananalysis point, similarity measurements with different dips usuallyrender a maximum similarity value in that case, where the dip of themeasurement coincides with the local dip in the data volume. In order tofind this dip dependent maximum with sufficient accuracy, the range ofpossible dips must be discretized sufficiently densely. This correspondsto a discretization of the solid angle.

[0113] A discretization of the solid angle is possible with a multitudeof point grids. A polar representation of such grids in the planerepresents the dip angle by the length of a radius vector, whereas theazimuth angle controls the direction of the radius vector. In thisrepresentation, a radial grid follows from the discretization of dip andazimuth angle in constant intervals. This uniform discretatization ofthe plane angles is advantageous for a statistical evaluation and forthe colour display of strike directions. However, it implies a stronglynon-uniform discretization of the solid angle. The quadratic andtriangular grid, on the contrary, render a very uniform discretizationof the solid angle. In order to observe a minimum accuracy of thediscretization, they require significantly less grid points than theradial grid. As a consequence, the radial grid requires a largercomputation effort.

[0114] According to the invention, a modified radial grid is proposedfor the discretization of the solid angle. The modified radial gridreduces the disadvantage of the non-uniform discretization of the solidangle in the radial grid, and the resulting computation effort, butgenerally preserves the discretization of the dip and azimuth angleswith constant steps, which is advantageous for further evaluations.

[0115] The discretization of the solid angle, and the calculation ofsimilarity values is carried out in one or two stages:

[0116] The first stage starts from a regular radial grid for thediscretization of the solid angle. The maximum dip angle of thediscretization is termed θ_(max). At a certain discrete dip angle θ_(i)with θ_(max)/2<θ_(i)≦θ_(max), a similarity value is calculated for everydiscretized azimuth. At dip angles θ_(i) with θ_(max)/4<θ_(i)≦θmax/2,the calculation is performed for every second azimuth, at dip anglesθ_(i) with θ_(max)/8<θ_(i)≦θ_(max)/4 for every forth azimuth, at dipangles θ_(i) with θ_(max)/16<θ_(i)≦θ_(max)/8 for every eighth azimuth,etc. The user of the method can define the factor, up to which theazimuth interval is coarsened. In our applications, a factor 4 hasproved to be effective. From the calculated similarity values, themaximum is determined, and the corresponding dip and azimuth angles areused as the first dip estimate. In an optional second stage, theaccuracy of the first estimate of the azimuth angle can be refined, ifthis azimuth estimation was performed with a coarsened interval. Thecoarsened azimuth intervals on both sides of the first estimate are thensubdivided by the original azimuth step, and corresponding similarityvalues are calculated as long as an increase of the similarity values isobserved. If the maximum of the similarity values from this second stageexceeds the maximum from the first stage, then the first dip estimate isreplaced by the dip and azimuth angles of this general maximum.

[0117] For the horizontal orientation of the weighted environment at thedip θ_(o)=0°, the present invention implies a multiple discretization ofthe orientation according to several azimuth angles, if the weightedenvironment exhibits a direction of preference. The reason for this isthe special definition of the azimuth angle, which determines both therotation of the environment in the horizontal plane, and the dipdirection, in case of an environment with a direction of preference. Inabsence of a direction of preference, a single discretization of thehorizontal orientation θ_(o)0° is sufficient.

[0118]FIG. 1 shows a modified radial grid in polar representation. Inthe first stage of the dip dependent calculation of similarity values,the modified grid uses 67 dips, which are marked as black grid points.Grid points, that are omitted in comparison to the original radial grid,are displayed in white.

[0119] In the second stage, the dip estimate of the first stage isoptionally refined to the accuracy of the original radial grid, as shownin FIG. 1. The assumed dip with maximum similarity of the first stage ismarked by a large black point. The dip angle of the maximum correspondsto a circle in this scheme. On this circle, the grid points of theoriginal radial scheme that have not been considered, are marked by sixlarge white points in the azimuth intervals on both sides of themaximum. Similarity values are calculated for these six dips, as long asan increase of the similarity values is observed. Thus, the number ofconsidered dips, or grid points, respectively, can increase to 73 in themodified radial scheme of FIG. 1.

[0120] In the present invention, this is prescribed for anon-symmetrical, weighted environment with a direction of preference.The reason for this is the deviating definition of the azimuth angle,which determines both the orientation of the environment in thehorizontal plane, and the dip direction, in case of an environment witha direction of preference.

[0121] Results

[0122] After the detection of t he dip dependent coherency maximum, thiscoherency maximum and the corresponding dip and azimuth angles arestored as results of the similarity measurements. They are thenavailable for a further digital evaluation, e.g. by pattern recognition.Furthermore, they are graphically displayed to the interpreter on themonitor screen, along cuts through the three-dimensional volume, oralong the processed time slices, or horizons, respectively. Suchdisplays are also printed on the printers that belong to the dataprocessing system.

APPLICATIONS EXAMPLE

[0123] For a three-dimensional volume of seismic poststack data, asimilarity volume was computed with help of the previously describedinvention. A weighted environment was defined by the Gaussian bellshaped exponential curve e^(−d) ² , where d is a scaled distanceaccording to equation (5b). A distance scaling according to a⁻¹=7.5 Δuwas used, with Δu being the average interval of the horizontaldiscretization in the seismic volume. Beyond a radius of 11.5 gridpoints around each analysis point, the weighting was set to the value 0.

[0124] The local similarity measurement in such an environment wasperformed with the weighted semblance according to equations (6b) and(10). The similarity measurements were executed in dipping environments.On principle, the dips were selected according to the scheme of FIG. 1.The spatial orientation is given by annotations N, denoting north, andE, denoting east. For each analysis position, the maximum similarityvalue that was thus determined, and the corresponding azimuth and dipangles were stored. These three parameters are displayed in FIG. 2a-cfor a horizontal plane, i.e., for a time slice of the seismic volume.

[0125] The value ranges are supplied with the corresponding gray scale.In FIG. 2a, the coherency is displayed in the value range from 0.0(black) to 1.0 (white) for a time slice from a seismic data volume. FIG.2b shows the local dip angle at maximum coherency for the coherencydisplay of FIG. 2a. The value range extends from −50° (black) to 0°(white) and +50° (black). FIG. 2c shows the derived local dip azimuthfor the coherency of FIG. 2a, and for the dip angle of FIG. 2b,respectively. The displayed value range extends from 0° (black) over 90°(white) to 180° (black).

1. A method for the determination of local similarity values forgeological units in the subsurface from a seismic 3-D dataset, whichconsists of a multitude of traces, which are formed in each case by asequence of data points, to which amplitude values, or seismicattributes derived from the amplitude values, are assigned, with thesteps: assigning an environment that is weighted at least horizontally,to a respective analysis point, whereby weight factors according to aweight function are assigned to the data points in the environment ofthe analysis point; calculating similarity values in each weightedenvironment, whereby the weight function is included in the similaritymeasure.
 2. A method for the determination of local dip dependentsimilarity values for geological units in the subsurface from a seismic3-D dataset, which consists of a multitude of traces, which are formedin each case by a sequence of data points, to which amplitude values, orseismic attributes derived from the amplitude values, are assigned, withthe steps: assigning a weighted environment to a respective analysispoint, whereby weight factors according to a weight function areassigned to the data points in the environment of the analysis point;calculating similarity values in each weighted environment, whereby theweight function is included in the similarity measure, and similarityvalues are determined for discrete spatial orientations, and thesimilarity values are represented at least as a function of a dip angleand a dip azimuth; determining, in each weighted environment of ananalysis point, the largest similarity value, which is assigned,together with the corresponding dip angle and dip azimuth, to therespective analysis point.
 3. The method according to claim 1, whereinthe weight function is centered at the analysis point, and exhibitsdecreasing values away from the analysis point.
 4. The method accordingto claim 1, wherein the weighting is embodied by window functions,especially as triangular, Hamming-, Hanning- or Daniell function.
 5. Themethod according to claim 1, wherein functions which continuously decaywith distance (d) up to infinify from the analysis point, especiallyd^(−c), e^(−cd), e^(−cdd) with c=arbitrary constant, are used for theweighting, whereby the functions have the value zero beyond a maximumdistance (d_(max)) from the analysis point.
 6. The method according toclaim 1, wherein the weighting is formed as the product of twofunctions, whereby one of the functions does not, or not completely,decay with the distance from the analysis point, and the other functionis a window function, which possesses the value 1 in the central area ofthe environment, comprises a decrease from 1 to 0 in the edge area, andadopts the value 0 outside the environment.
 7. The method according toclaim 1, wherein the weighted environment is formed with rotationalsymmetry to the analysis point at least in the horizontal plane.
 8. Themethod according to claim 2, with one of the claims 3 to 6, wherein theweighted environment, which is not formed with rotational symmetry tothe analysis point in the horizontal plane, and which exhibits adirection of preference, is rotated with the respective dip azimuth, forthe determination of similarity values.
 9. The method according to claim1, wherein the interval between neighboring analysis points in the 3-Ddataset is selected as multiples, and especially as integer multiples ofthe data point intervals, preferably between one and several tens,especially preferred between three and ten data point intervals.
 10. Themethod according to claim 2, with one of the claims 3 to 9, wherein thedip and azimuth angles in a radial grid are each discretized in constantangle steps, and that, in the calculation of the similarity values withdependence on the discrete dip angle (θ_(i)), starting from apre-defined maximum dip angle (θ_(max)), the resulting discrete dips areconsidered as follows: at dip angles θ_(i) with θ_(max)/2<θ_(i)≦θ_(max),a similarity value is calculated for every discretized azimuth, at dipangles θ_(i) with θ_(max)/4<θ_(i)≦θ_(max)/2, the calculation isperformed for every second azimuth, at dip angles θ_(i) withθ_(max)/8<θ_(i)≦θ_(max)/4 for every forth azimuth, at dip angles θ_(i)with θ_(max)/16<θ_(i)≦θ_(max)/8 for every eighth azimuth, etc., with apreferred termination at the fourth subdivision.
 11. The methodaccording to claim 10, wherein for a dip estimate in a coarsened azimuthinterval, similarity values are calculated at the original angle stepson both sides of the first estimate, as long as an increase of thesimilarity values is observed, whereby the first dip estimate isreplaced by the dip and azimuth angles of the this general similaritymaximum, if the maximum of the similarity values from this second stageexceeds the maximum from the first stage.
 12. The method according toclaim 1, wherein the similarity is computed along arbitrarily orientedplanes as the square of the representative attribute (amplitude),divided by the average of the squared attribute (amplitude).
 13. Themethod according to claim 12, wherein the representative attribute isthe median of the attributes to be analyzed.
 14. The method according toclaim 13, wherein horizontal data slice at a time or depthz_({overscore (k)}) with a total of N amplitude valuess(x_(i),y_(j),z_({overscore (k)})), as well as N corresponding weightfactors g_(IJK)(x_(i),y_(j),z_({overscore (k)})) is considered, wherebythe amplitude values are sorted in increasing order, and newly indexed,as: s₁≦s₂≦s₃≦ . . . ≦s_(N−1)≦s_(N), this sorting and indexing istransferred to the values of the weight function g_(IJK) belonging tothe seismic amplitudes, as: g₁, g₂, g₃, . . . g_(N−1), g_(N), andpartial sums of the weight factors are derived G₁≦G₂≦G₃≦ . . .≦G_(N−1)≦G_(N), with G_(n)=Σg_(i)(i=1 to n), wherefrom the weightedmedian can be calculated as${m_{IJK}( z_{\overset{\_}{k}} )} = \begin{Bmatrix}{s_{n},} & {{{{if}\quad G_{n - 1}} < \frac{G_{N}}{2}},{G_{n} > \frac{G_{N}}{2}}} \\{\frac{s_{n} + s_{m}}{2},} & {{{{if}\quad G_{n - 1}} < \frac{G_{N}}{2}},{G_{n} = {\frac{G_{N}}{2} = G_{m - 1}}},{G_{m} > \frac{G_{N}}{2}}} \\\quad & {{{with}\quad n} < m}\end{Bmatrix}$

and from this, the weighted median similarity follow for a horizontaldata slice at an arbitrary time or depth z_({overscore (k)}):${M_{\overset{\_}{k}}( {x_{I},y_{J},z_{K}} )} = \frac{\{ {\sum\limits_{i,j}{g_{IJK}( {x_{i},y_{j},z_{\overset{\_}{k}}} )}} \} \{ {m_{IJK}( z_{\overset{\_}{k}} )} \}^{2}}{\sum\limits_{i,j}{{g_{IJK}( {x_{i},y_{j},z_{\overset{\_}{k}}} )}\{ {s( {x_{i},y_{j},z_{\overset{\_}{k}}} )} \}^{2}}}$

wherefrom, by summation over the entire time or depth range of theweighted environment, the weighted median coherency follows as asimilarity measure, as:${M( {x_{I},y_{J},z_{K}} )} = \frac{\sum\limits_{k}\lbrack {\{ {\sum\limits_{i,j}{g_{IJK}( {x_{i},y_{j},z_{k}} )}} \} \{ {m_{IJK}( z_{k} )} \}^{2}} \rbrack}{\sum\limits_{k}\lbrack {\sum\limits_{i,j}{{g_{IJK}( {x_{i},y_{j},z_{k}} )}\{ {s( {x_{i},y_{j},z_{k}} )} \}^{2}}} \rbrack}$


15. The method according to claim 1, wherein the similarity iscalculated as semblance modified to consider the weighted environmentg_(IJK)(x_(i),y_(j),z_(k)), for an analysis position (x_(I),y_(J),z_(K))as${{S_{\overset{\_}{k}}( {x_{I},y_{J},z_{K}} )} = \frac{\{ {\sum\limits_{i,j}{g_{IJK}( {x_{i},y_{j},z_{\overset{\_}{k}}} )}} \}^{- 1}\{ {\sum\limits_{i,j}{{g_{IJK}( {x_{i},y_{j},z_{\overset{\_}{k}}} )}{s( {x_{i},y_{j},z_{\overset{\_}{k}}} )}}} \}^{2}}{\sum\limits_{i,j}{{g_{IJK}( {x_{i},y_{j},z_{\overset{\_}{k}}} )}\{ {s( {x_{i},{y_{j}z_{\overset{\_}{k}}}} )} \}^{2}}}},$

wherefrom, by summation over the entire time or depth range of theweighted environment, the similarity measure of the weighted coherencyfollows:${S( {x_{I},y_{J},z_{K}} )} = \frac{\sum\limits_{k}\lbrack {\{ {\sum\limits_{i,j}{g_{IJK}( {x_{i},y_{j},z_{k}} )}} \}^{- 1}\{ {\sum\limits_{i,j}{{g_{IJK}( {x_{i},y_{j},z_{k}} )}{s( {x_{i},y_{j},z_{k}} )}}} \}^{2}} \rbrack}{\sum\limits_{k}\lbrack {\sum\limits_{i,j}{{g_{IJK}( {x_{i},y_{j},z_{k}} )}\{ {s( {x_{i},{y_{j}z_{k}}} )} \}^{2}}} \rbrack}$


16. The method according to claim 2, wherein the weight function iscentered at the analysis point, and exhibits decreasing values away fromthe analysis point.
 17. The method according to claim 2, wherein theweighting is embodied by window functions, especially as triangular,Hamming-, Hanning- or Daniell function.
 18. The method according toclaim 2, wherein functions which continuously decay with distance (d) upto infinify from the analysis point, especially d^(−c), e^(−cd),e^(−cdd) with c=arbitrary constant, are used for the weighting, wherebythe functions have the value zero beyond a maximum distance (d_(max))from the analysis point.
 19. The method according to claim 2, whereinthe weighting is formed as the product of two functions, whereby one ofthe functions does not, or not completely, decay with the distance fromthe analysis point, and the other function is a window function, whichpossesses the value 1 in the central area of the environment, comprisesa decrease from 1 to 0 in the edge area, and adopts the value 0 outsidethe environment.
 20. The method according to claim 2, wherein theweighted environment is formed with rotational symmetry to the analysispoint at least in the horizontal plane.